clear all
set more off, perm
set mem 10000000
set matsize 10000

********************************************************* 
*** DD regressions: first-stage nighttime brightness ****
********************************************************* 

** Set file paths
do "$path_code/paths.do"

** Set graph scheme
cd "$path/code/analyze"
set scheme fb, perm

********************************************************************************
********************************************************************************

** 1. DD pooled, and split by population bin (1998-2013)
{
use "${panel}/panel_dataset_full.dta", clear
keep pca01_id st_code dt_code state district stdt vplan4 vdpr4 ///
	lights_mean???? lights_max???? lights_sta_max???? lights_pct_max????
reshape long lights_mean lights_max lights_sta_max lights_pct_max, i(pca01_id) j(year)
la var year "Year"
la var lights_mean "Mean village brightness"
la var lights_max "Max village brightness"
la var lights_sta_max "Max village brightness (STABLE)"
la var lights_pct_max "Max village brightness (PCT)"
compress
save "${panel}/panel_dataset_lights_allyears.dta", replace


use "${panel}/panel_dataset_dd.dta", clear

	// drop 10th-Plan states with bad shapefiles
drop if inlist(state,"UTTARAKHAND","UTTAR PRADESH","HIMACHAL PRADESH","ASSAM","JAMMU & KASHMIR")
	
	// drop states with missing shapefiles
drop if inlist(state,"ARUNACHAL PRADESH","MEGHALAYA","MIZORAM","NAGALAND","SIKKIM")

	// drop 11th-Plan state with bad shapefile
drop if inlist(state,"PUNJAB")
	
	// drop if lights are missing
drop if lights_max==.
tab vplan4 sample_1011, missing

drop pct_06-vd_ame_d_spcl
drop vd_com_d_newspap-vd_ame_d_cv_hall
drop work_main_ag_p-pct_irr_twell pct_irr_irr_twell
drop p_300_X-band300_p_200
drop bin25_p-within_border_100k
drop vplan vplan2 vplan3 vdpr vdpr2 vdpr3
drop treat_x_post dd_* stateyear
drop lights_mean_hat lights_max_hat
drop no_hh tot_p tot_m tot_f vd01_id conc_id v_id_hpca11 vd11_id c_code01 v_id_master ///
	h?_id h?v_id f_area totp_h? count_h? maxp_h? grad_max grad_mean
merge 1:1 pca01_id st_code dt_code vplan4 vdpr4 year using ///
	"${panel}/panel_dataset_lights_allyears.dta"
egen max_merge = max(_merge), by(pca01_id)
keep if max_merge==3
tab year _merge
assert _merge==3 if year==2001 | year==2011
drop _merge max_merge

order pca01_id year state st_code district dt_code stdt vplan4 vdpr4 lights* vd_pwr*
preserve
keep if inlist(year,2001,2011)
keep pca01_id bk_code-popbin_2700_max
duplicates drop
unique pca01_id
assert r(unique)==r(N)
tempfile vars_constant
save `vars_constant'
restore
drop bk_code-popbin_2700_max
merge m:1 pca01_id using `vars_constant'
assert _merge==3
drop _merge
sort pca01_id year


preserve
use "$rggvy/rggvy_district_progress_X_XI_processed.dta", clear
replace award_date = sanction_date if award_date<sanction_date
replace award_date = max(award_date,17553) if plan==11
collapse (min) min_award_date=award_date (max) max_award_date=award_date, ///
	by(st_code dt_code) fast
gen med_award_date = round((min_award_date+max_award_date)/2,1)
format %td med_award_date
tempfile award_dates
save `award_dates'
restore
merge m:1 st_code dt_code using `award_dates', keep(1 3)
assert _merge==3 if vplan4!=.
assert _merge==1 if vplan4==.
drop _merge

preserve
use "$panel/panel_dataset_dd_nss.dta", clear
keep st_code dt_code exp05*
duplicates drop
tempfile nss
save `nss'
restore
merge m:1 st_code dt_code using `nss', nogen keep(1 3)


	// DD 1: 10th plan * post-2005
gen dd_1_treat_x_post = vplan4<11 & year>2005
la var dd_1_treat_x_post "DD 1: post-2005 (both plans)"	

	// DD 2: 10th plan * post-award-year
gen award_year = year(med_award_date)
egen stdt_tag = tag(st_code dt_code)
tab award_year vplan4 if stdt_tag, missing
gen dd_2_treat_x_post = vplan4<11 & year>award_year
la var dd_2_treat_x_post "DD 2: post-award-year (10th plan)"	

	// DD 3: post-award-year
gen dd_3_treat_x_post = year>award_year
la var dd_3_treat_x_post "DD 3: post-award-year (both plans)"	
	
	// Create 500-person population bins (switching from 300 to 500 to align with NSS population medians)
drop popbin*	
gen popbin = 0
local bin = 500
local bmax = 3000
local bmax2 = 2900	
forvalues i = 0(`bin')`bmax2' {
  local iplus1 = `i' + `bin'
  gen popbin_`i'_`iplus1' = tot_p01>`i' & tot_p01<=`iplus1'
  la var popbin_`i'_`iplus1' "2001 population bin: (`i',`iplus1']"
  replace popbin = `i' if popbin_`i'_`iplus1'==1
}
gen popbin_`bmax'_max = tot_p01>`bmax'
replace popbin = `bmax' if popbin_`bmax'_max==1
la var popbin "Lower end of 2001 population bin"

	// DD population bins
foreach j in 1 2 3 {
	forvalues i = 0(`bin')`bmax2' {
	  local iplus1 = `i' + `bin'
	  gen dd_`j'_x_popbin_`i'_`iplus1' = dd_`j'_treat_x_post * popbin_`i'_`iplus1'
	  la var dd_`j'_x_popbin_`i'_`iplus1' "DD `j' treat x popbin: `i'-`iplus1'"
	}
	gen dd_`j'_x_popbin_`bmax'_max = dd_`j'_treat_x_post * popbin_`bmax'_max
	la var dd_`j'_x_popbin_`bmax'_max "DD `j' treat x popbin `bmax' to infinity"
}


	// REGRESSION LOOP
cap drop yvar-rmse
cap drop row_id
gen yvar = ""
gen dd_model = ""
gen dd_var = .
gen fes = ""
gen ifs = ""
gen beta = .
gen se = .
gen tscore = .
gen pvalue = .
gen ci95_lo = .
gen ci95_hi = .
forvalues p = 0(`bin')`bmax' {
	local p1 = `p'+`bin'
	if `p'==`bmax' {
		local p1 = "max"
	}
	gen beta_`p'_`p1' = .
	gen se_`p'_`p1' = .
	gen tscore_`p'_`p1' = .
	gen pvalue_`p'_`p1' = .
	gen ci95_lo_`p'_`p1' = .
	gen ci95_hi_`p'_`p1' = .
}
gen ymean = .
gen nobs = .
gen nclust = .
gen rmse = .
gen row_id = .

local row = 1

foreach dd_var in 1 2 3 {
	
	foreach dd_sample in 1 2 {
		
		if `dd_sample'==1 {
			local ifs = ""
		}
		else if `dd_sample'==2 {
			local ifs = " if sample_1011==1"
		}
		
		foreach yvar of varlist lights_max {
	
			foreach FEs in 1 2 3 {
				
				if `FEs'==1 {
					local fes = "pca01_id year"
				}
				else if `FEs'==2 {
					local fes = "pca01_id year exp05_st_4ile#c.year exp05_ntl_10ile#c.year"
				}
				else if `FEs'==3 {
					local fes = "pca01_id year exp05_st_4ile#c.year exp05_ntl_10ile#c.year c.year#st_code"
				}
				
				*pooled DD
				local dd_model = "pooled"
				reghdfe `yvar' dd_`dd_var'_treat_x_post `ifs', absorb(`fes') vce(cluster stdt)
				qui replace yvar = "`yvar'"	if _n==`row'
				qui replace dd_model = "`dd_model'" if _n==`row'
				qui replace dd_var = `dd_var'	if _n==`row'
				qui replace fes = "`fes'" if _n==`row'
				qui replace ifs = "`ifs'" if _n==`row'
				qui replace beta = _b[dd_`dd_var'_treat_x_post] if _n==`row'
				qui replace se = _se[dd_`dd_var'_treat_x_post]	if _n==`row'
				qui replace tscore = beta/se in `row'
				qui replace pvalue = 2*ttail(e(df_r),abs(_b[dd_`dd_var'_treat_x_post]/_se[dd_`dd_var'_treat_x_post])) in `row'
				qui replace ci95_lo = beta - se*invttail(e(df_r),0.025) in `row'
				qui replace ci95_hi = beta + se*invttail(e(df_r),0.025) in `row'
				qui sum `yvar' if e(sample)	
				qui replace ymean = r(mean)	if _n==`row'
				qui replace nobs = e(N)	if _n==`row'
				qui replace nclust = e(N_clust)	if _n==`row'
				qui replace rmse = e(rmse) if _n==`row'
				qui replace row_id	= `row'	if _n==`row'
				qui local row = `row' + 1
				
				* DD w/ pop bins
				local dd_model = "pop bins"
				local fes = subinstr("`fes'"," year ", " year#popbin ",1)
				reghdfe `yvar' dd_`dd_var'_x_popbin_*  `ifs', absorb(`fes') vce(cluster stdt)
				qui replace yvar = "`yvar'"	if _n==`row'
				qui replace dd_model = "`dd_model'" if _n==`row'
				qui replace dd_var = `dd_var'	if _n==`row'
				qui replace fes = "`fes'" if _n==`row'
				qui replace ifs = "`ifs'" if _n==`row'
				forvalues p = 0(`bin')`bmax' {	
					local p1 = `p'+`bin'
					if `p'==`bmax' {
						local p1 = "max"
					}
					qui replace beta_`p'_`p1' = _b[dd_`dd_var'_x_popbin_`p'_`p1'] if _n==`row'
					qui replace se_`p'_`p1' = _se[dd_`dd_var'_x_popbin_`p'_`p1'] if _n==`row'
					qui replace tscore_`p'_`p1' = beta_`p'_`p1'/se_`p'_`p1' in `row'
					qui replace pvalue_`p'_`p1' = 2*ttail(e(df_r),abs(_b[dd_`dd_var'_x_popbin_`p'_`p1']/_se[dd_`dd_var'_x_popbin_`p'_`p1'])) in `row'
					qui replace ci95_lo_`p'_`p1' = beta_`p'_`p1' - se_`p'_`p1'*invttail(e(df_r),0.025) in `row'
					qui replace ci95_hi_`p'_`p1' = beta_`p'_`p1' + se_`p'_`p1'*invttail(e(df_r),0.025) in `row'
				}
				qui sum `yvar' if e(sample)							
				qui replace ymean = r(mean)	if _n==`row'
				qui replace nobs = e(N)	if _n==`row'
				qui replace nclust = e(N_clust)	if _n==`row'
				qui replace rmse = e(rmse)	if _n==`row'
				qui replace row_id	= `row'	if _n==`row'
				qui local row = `row' + 1
			}
		}
	}
}

	// save results
keep yvar-row_id
keep if yvar!=""
compress
save "$results/dd_results_lights_allyears.dta", replace

}

********************************************************************************
********************************************************************************

** 2. Event-study DD
{
use "${panel}/panel_dataset_dd.dta", clear
keep if sample==1 & sample_1011==1 & corr_state==1
drop pct_06-vd_ame_d_spcl
drop vd_com_d_newspap-vd_ame_d_cv_hall
drop work_main_ag_p-pct_irr_twell pct_irr_irr_twell
drop p_300_X-band300_p_200
drop bin25_p-within_border_100k
drop vplan vplan2 vplan3 vdpr vdpr2 vdpr3
drop treat_x_post dd_* stateyear
drop lights_mean_hat lights_max_hat
drop no_hh tot_p tot_m tot_f vd01_id conc_id v_id_hpca11 vd11_id c_code01 v_id_master ///
	h?_id h?v_id f_area totp_h? count_h? maxp_h? grad_max grad_mean
merge 1:1 pca01_id st_code dt_code vplan4 vdpr4 year using ///
	"${panel}/panel_dataset_lights_allyears.dta"
egen max_merge = max(_merge), by(pca01_id)
keep if max_merge==3
tab year _merge
assert _merge==3 if year==2001 | year==2011
drop _merge max_merge

order pca01_id year state st_code district dt_code stdt vplan4 vdpr4 lights* vd_pwr*
preserve
keep if inlist(year,2001,2011)
keep pca01_id bk_code-popbin_2700_max
duplicates drop
unique pca01_id
assert r(unique)==r(N)
tempfile vars_constant
save `vars_constant'
restore
drop bk_code-popbin_2700_max
merge m:1 pca01_id using `vars_constant'
assert _merge==3
drop _merge
sort pca01_id year


preserve
use "$rggvy/rggvy_district_progress_X_XI_processed.dta", clear
replace award_date = sanction_date if award_date<sanction_date
replace award_date = max(award_date,17553) if plan==11
collapse (min) min_award_date=award_date (max) max_award_date=award_date, ///
	by(st_code dt_code) fast
gen med_award_date = round((min_award_date+max_award_date)/2,1)
format %td med_award_date
tempfile award_dates
save `award_dates'
restore
merge m:1 st_code dt_code using `award_dates', keep(1 3)
assert _merge==3
drop _merge
egen stdtbk = group(stdt bk_code)
*egen lights_max_dist = mean(lights_max), by(stdt year)

	// DD 1: 10th plan * post-2005 (year zero)
gen dd_1_treat_x_post = vplan4==10 & year==2005
la var dd_1_treat_x_post "DD 1: post-2005 (both plans)"	

	// DD 2: 10th plan * post-award-year (year zero)
gen award_year = year(med_award_date)
assert award_year!=.	
egen stdt_tag = tag(st_code dt_code)
tab award_year vplan4 if stdt_tag
gen dd_2_treat_x_post = vplan4==10 & year==award_year
la var dd_2_treat_x_post "DD 2: post-award-year (10th plan)"	

	// DD 3: post-award-year (year zero)
gen dd_3_treat_x_post = year==award_year
la var dd_3_treat_x_post "DD 3: post-award-year (both plans)"	
	
	// Set time series
tsset pca01_id year	
	
	// Populate leads and lags
cap drop dd_?_pre? dd_?_t0 dd_?_post?
foreach v in 3 2 {
	gen dd_`v'_pre13 = F13.dd_`v'_treat_x_post
	gen dd_`v'_pre12 = F12.dd_`v'_treat_x_post
	gen dd_`v'_pre11 = F11.dd_`v'_treat_x_post
	gen dd_`v'_pre10 = F10.dd_`v'_treat_x_post
	gen dd_`v'_pre9 = F9.dd_`v'_treat_x_post
	gen dd_`v'_pre8 = F8.dd_`v'_treat_x_post
	gen dd_`v'_pre7 = F7.dd_`v'_treat_x_post
	gen dd_`v'_pre6 = F6.dd_`v'_treat_x_post
	gen dd_`v'_pre5 = F5.dd_`v'_treat_x_post
	gen dd_`v'_pre4 = F4.dd_`v'_treat_x_post
	gen dd_`v'_pre3 = F3.dd_`v'_treat_x_post
	gen dd_`v'_pre2 = F2.dd_`v'_treat_x_post
	gen dd_`v'_pre1 = F1.dd_`v'_treat_x_post
	gen dd_`v'_t0 = dd_`v'_treat_x_post
	gen dd_`v'_post1 = L1.dd_`v'_treat_x_post
	gen dd_`v'_post2 = L2.dd_`v'_treat_x_post
	gen dd_`v'_post3 = L3.dd_`v'_treat_x_post
	gen dd_`v'_post4 = L4.dd_`v'_treat_x_post
	gen dd_`v'_post5 = L5.dd_`v'_treat_x_post
	gen dd_`v'_post6 = L6.dd_`v'_treat_x_post
	gen dd_`v'_post7 = L7.dd_`v'_treat_x_post
	gen dd_`v'_post8 = L8.dd_`v'_treat_x_post
	gen dd_`v'_post9 = L9.dd_`v'_treat_x_post
}	
*br pca01_id year vplan4 award_year dd_1_*	
foreach v in  3 2 {
	egen temp = rowmax(dd_`v'_post5 dd_`v'_post6 dd_`v'_post7 dd_`v'_post8 dd_`v'_post9)
	replace dd_`v'_post5 = temp
	drop dd_`v'_post6 dd_`v'_post7 dd_`v'_post8 dd_`v'_post9 temp
}
foreach v in 3 2 {
	egen temp = rowmax(dd_`v'_pre13 dd_`v'_pre12 dd_`v'_pre11 dd_`v'_pre10 dd_`v'_pre9 dd_`v'_pre8 dd_`v'_pre7 dd_`v'_pre6 dd_`v'_pre5)
	replace dd_`v'_pre5 = temp
	drop dd_`v'_pre13 dd_`v'_pre12 dd_`v'_pre11 dd_`v'_pre10 dd_`v'_pre9 dd_`v'_pre8 dd_`v'_pre7 dd_`v'_pre6 temp
}

foreach v in 3 2 {
	foreach v2 of varlist dd_`v'_pre* dd_`v'_t0 dd_`v'_post* {
		replace `v2' = 0 if `v2'==.
	}
}

	// Event study using funding from both Plans
reghdfe lights_max dd_3_pre5-dd_3_pre2 dd_3_t0 dd_3_post* , absorb(pca01_id year st_code#year pca01_id#c.year) vce(cluster stdtbk)

cap drop t beta lci uci
gen t = _n-6
replace t= . if t>5	
gen beta = .
gen lci = .
gen uci = .
replace beta = _b[dd_3_pre5] if t==-5	
replace beta = _b[dd_3_pre4] if t==-4	
replace beta = _b[dd_3_pre3] if t==-3	
replace beta = _b[dd_3_pre2] if t==-2
cap replace beta = _b[dd_3_pre1] if t==-1
if _rc {
	replace beta = 0 if t==-1
}
cap replace beta = _b[dd_3_t0] if t==0
if _rc {
	replace beta = 0 if t==0
}
replace beta = _b[dd_3_post1] if t==1
replace beta = _b[dd_3_post2] if t==2
replace beta = _b[dd_3_post3] if t==3
replace beta = _b[dd_3_post4] if t==4
replace beta = _b[dd_3_post5] if t==5

replace lci = _b[dd_3_pre5]-1.96*_se[dd_3_pre5] if t==-5	
replace lci = _b[dd_3_pre4]-1.96*_se[dd_3_pre4] if t==-4	
replace lci = _b[dd_3_pre3]-1.96*_se[dd_3_pre3] if t==-3	
replace lci = _b[dd_3_pre2]-1.96*_se[dd_3_pre2]  if t==-2
cap replace lci = _b[dd_3_pre1]-1.96*_se[dd_3_pre1]  if t==-1
cap replace lci = _b[dd_3_t0]-1.96*_se[dd_3_t0]  if t==0
replace lci = _b[dd_3_post1]-1.96*_se[dd_3_post1]  if t==1
replace lci = _b[dd_3_post2]-1.96*_se[dd_3_post2]  if t==2
replace lci = _b[dd_3_post3]-1.96*_se[dd_3_post3]  if t==3
replace lci = _b[dd_3_post4]-1.96*_se[dd_3_post4]  if t==4
replace lci = _b[dd_3_post5]-1.96*_se[dd_3_post5]  if t==5

replace uci = _b[dd_3_pre5]+1.96*_se[dd_3_pre5] if t==-5	
replace uci = _b[dd_3_pre4]+1.96*_se[dd_3_pre4] if t==-4	
replace uci = _b[dd_3_pre3]+1.96*_se[dd_3_pre3] if t==-3	
replace uci = _b[dd_3_pre2]+1.96*_se[dd_3_pre2]  if t==-2
cap replace uci = _b[dd_3_pre1]+1.96*_se[dd_3_pre1]  if t==-1
cap replace uci = _b[dd_3_t0]+1.96*_se[dd_3_t0]  if t==0
replace uci = _b[dd_3_post1]+1.96*_se[dd_3_post1]  if t==1
replace uci = _b[dd_3_post2]+1.96*_se[dd_3_post2]  if t==2
replace uci = _b[dd_3_post3]+1.96*_se[dd_3_post3]  if t==3
replace uci = _b[dd_3_post4]+1.96*_se[dd_3_post4]  if t==4
replace uci = _b[dd_3_post5]+1.96*_se[dd_3_post5]  if t==5
	
twoway (rspike uci lci t, lstyle(ci) lw(medium) lcolor(gs10)) ///
	(scatter beta t, msize(medsmall) mc(gs4) msymbol(O)) ///
	, ///
	yline(0, lcolor(black) lw(thin)) /// 
	ytitle("Event-study  coefficient (95% CI)", size(medlarge)) xtitle("Event time", size(medlarge)) ///
	title("Treatment: when any district received funds", color(black) size(large)) /// 
	graphregion(lcolor(white)) graphregion(color(white)) plotregion(fcolor(white)) ///
	ylabel(/*-0.03 -0.02 -0.01 0 0.01*/,nogrid angle(0) labsize(medlarge)) /*yscale(range(-0.015 0.015))*/ ///
	xlabel(-5 "-5+" -4 -3 -2 -1 0 1 2 3 4 5 "5+", labsize(medlarge)) ///
	 legend(off)

graph export "$results/dd figs/event_study_lights_bothplans.eps", replace 


	// Event study using funding from 10th Plan only
reghdfe lights_max dd_2_pre5-dd_2_pre2 dd_2_t0 dd_2_post* , absorb(pca01_id year st_code#year pca01_id#c.year) vce(cluster stdtbk)

cap drop t beta lci uci
gen t = _n-6
replace t= . if t>5	
gen beta = .
gen lci = .
gen uci = .
replace beta = _b[dd_2_pre5] if t==-5	
replace beta = _b[dd_2_pre4] if t==-4	
replace beta = _b[dd_2_pre3] if t==-3	
replace beta = _b[dd_2_pre2] if t==-2
cap replace beta = _b[dd_2_pre1] if t==-1
if _rc {
	replace beta = 0 if t==-1
}
cap replace beta = _b[dd_2_t0] if t==0
if _rc {
	replace beta = 0 if t==0
}
replace beta = _b[dd_2_post1] if t==1
replace beta = _b[dd_2_post2] if t==2
replace beta = _b[dd_2_post3] if t==3
replace beta = _b[dd_2_post4] if t==4
replace beta = _b[dd_2_post5] if t==5

replace lci = _b[dd_2_pre5]-1.96*_se[dd_2_pre5] if t==-5	
replace lci = _b[dd_2_pre4]-1.96*_se[dd_2_pre4] if t==-4	
replace lci = _b[dd_2_pre3]-1.96*_se[dd_2_pre3] if t==-3	
replace lci = _b[dd_2_pre2]-1.96*_se[dd_2_pre2]  if t==-2
cap replace lci = _b[dd_2_pre1]-1.96*_se[dd_2_pre1]  if t==-1
cap replace lci = _b[dd_2_t0]-1.96*_se[dd_2_t0]  if t==0
replace lci = _b[dd_2_post1]-1.96*_se[dd_2_post1]  if t==1
replace lci = _b[dd_2_post2]-1.96*_se[dd_2_post2]  if t==2
replace lci = _b[dd_2_post3]-1.96*_se[dd_2_post3]  if t==3
replace lci = _b[dd_2_post4]-1.96*_se[dd_2_post4]  if t==4
replace lci = _b[dd_2_post5]-1.96*_se[dd_2_post5]  if t==5

replace uci = _b[dd_2_pre5]+1.96*_se[dd_2_pre5] if t==-5	
replace uci = _b[dd_2_pre4]+1.96*_se[dd_2_pre4] if t==-4	
replace uci = _b[dd_2_pre3]+1.96*_se[dd_2_pre3] if t==-3	
replace uci = _b[dd_2_pre2]+1.96*_se[dd_2_pre2]  if t==-2
cap replace uci = _b[dd_2_pre1]+1.96*_se[dd_2_pre1]  if t==-1
cap replace uci = _b[dd_2_t0]+1.96*_se[dd_2_t0]  if t==0
replace uci = _b[dd_2_post1]+1.96*_se[dd_2_post1]  if t==1
replace uci = _b[dd_2_post2]+1.96*_se[dd_2_post2]  if t==2
replace uci = _b[dd_2_post3]+1.96*_se[dd_2_post3]  if t==3
replace uci = _b[dd_2_post4]+1.96*_se[dd_2_post4]  if t==4
replace uci = _b[dd_2_post5]+1.96*_se[dd_2_post5]  if t==5
	
twoway (rspike uci lci t, lstyle(ci) lw(medium) lcolor(gs10)) ///
	(scatter beta t, msize(medsmall) mc(gs4) msymbol(O)) ///
	, ///
	yline(0, lcolor(black) lw(thin)) /// 
	ytitle("Event-study  coefficient (95% CI)", size(medlarge)) xtitle("Event time", size(medlarge)) ///
	title("Treatment: when 1st-wave districts received funds", color(black) size(large)) /// 
	graphregion(lcolor(white)) graphregion(color(white)) plotregion(fcolor(white)) ///
	ylabel(/*-0.03 -0.02 -0.01 0 0.01*/,nogrid angle(0) labsize(medlarge))  yscale(titlegap(*-30)) ///
	xlabel(-5 "-5+" -4 -3 -2 -1 0 1 2 3 4 5 "5+", labsize(medlarge)) ///
	 legend(off)

graph export "$results/dd figs/event_study_lights_10thplan.eps", replace 

}

********************************************************************************
********************************************************************************

** 3. Difference-in-Discontinuities
{
use "${panel}/panel_dataset_dd.dta", clear
keep if vplan4<11 & corr_state==1 & sample==1 & sample_h==1 & sing_h==1 & pop_non_zero==1 
drop pct_06-vd_ame_d_spcl
drop vd_com_d_newspap-vd_ame_d_cv_hall
drop work_main_ag_p-pct_irr_twell pct_irr_irr_twell
drop p_300_X-band300_p_200
drop bin25_p-within_border_100k
drop vplan vplan2 vplan3 vdpr vdpr2 vdpr3
drop treat_x_post dd_* stateyear
drop lights_mean_hat lights_max_hat
drop no_hh tot_p tot_m tot_f vd01_id conc_id v_id_hpca11 vd11_id c_code01 v_id_master ///
	h?_id h?v_id f_area totp_h? count_h? maxp_h? grad_max grad_mean
merge 1:1 pca01_id st_code dt_code vplan4 vdpr4 year using ///
	"${panel}/panel_dataset_lights_allyears.dta"
egen max_merge = max(_merge), by(pca01_id)
keep if max_merge==3
tab year _merge
assert _merge==3 if year==2001 | year==2011
drop _merge max_merge

order pca01_id year state st_code district dt_code stdt vplan4 vdpr4 lights* vd_pwr*
preserve
keep if inlist(year,2001,2011)
keep pca01_id bk_code-popbin_2700_max
duplicates drop
unique pca01_id
assert r(unique)==r(N)
tempfile vars_constant
save `vars_constant'
restore
drop bk_code-popbin_2700_max
merge m:1 pca01_id using `vars_constant'
assert _merge==3
drop _merge
sort pca01_id year

preserve
use "${panel}/panel_dataset_full.dta", clear
keep pca01_id lights_max????_hat
rename lights_max????_hat lights_max_hat????
reshape long lights_max_hat, i(pca01_id) j(year)
tempfile lights_max_hat
save `lights_max_hat'
restore
merge 1:1 pca01_id year using `lights_max_hat', keep(1 3)
drop if tot_p01>1000

egen temp1 = max(lights_max*(year==2001)), by(pca01_id)
egen temp2 = max(lights_max*(year==2011)), by(pca01_id)
gen lights_diff = abs(temp2-temp1)
assert lights_diff!=.

egen stdtbk = group(stdt bk_code)
gen X = (tot_p01-300)
assert X!=.
gen D = X>=0
gen DX = D*X

gen aw150 = abs(150-abs(300-tot_p01))
gen aw100 = abs(100-abs(300-tot_p01))


 // District FEs, 100 bw
reghdfe lights_max ///
	 1.D#1998.year 1.D#1999.year 1.D#2000.year 1.D#2001.year 1.D#2002.year 1.D#2003.year 1.D#2004.year /// 1.D#2005.year ///
	 1.D#2006.year 1.D#2007.year 1.D#2008.year 1.D#2009.year 1.D#2010.year 1.D#2011.year 1.D#2012.year 1.D#2013.year ///
	 c.X#i.year c.DX#i.year if inrange(tot_p01,200,400) & pop_mismatch20==0 & lights_diff<20 [aw=aw100], a(pca01_id year#stdt) vce(cluster stdtbk) 
	 
cap drop y beta lci uci
gen y = _n+1997
replace y = . if y>2013
gen beta = .
gen lci = .
gen uci = .

forvalues y = 1998/2013 {
	cap replace beta = _b[1.D#`y'.year] if y==`y'
	cap replace lci =  _b[1.D#`y'.year] - 1.96*_se[1.D#`y'.year] if y==`y'
	cap replace uci =  _b[1.D#`y'.year] + 1.96*_se[1.D#`y'.year] if y==`y'
}
cap replace beta = 0 if y==2005

twoway (rspike uci lci y, lstyle(ci) lw(medium) lcolor(gs10)) ///
	(scatter beta y, msize(medsmall) mc(midblue) msymbol(O)) ///
	, ///
	yline(0, lcolor(black) lw(thin)) /// 
	ytitle("Year-specific RD coefficient (95% CI)", size(medlarge)) xtitle("Nighttime brightness year", size(medlarge)) ///
	title("District-year FEs, 100-person bandwidth", color(black) size(large)) /// 
	graphregion(lcolor(white)) graphregion(color(white)) plotregion(fcolor(white)) ///
	ylabel(-0.1 0 0.1 0.2 0.3 0.4,nogrid angle(0) labsize(medlarge)) yscale(range(-0.1 0.43)) ///
	xlabel(1998(1)2013, labsize(medlarge) angle(45)) ///
	 legend(off) //aspect(0.4)

graph export "$results/dd figs/diff_in_disc_lights_dt_100.pdf", replace 
	 

 // District FEs, 150 bw
reghdfe lights_max ///
	 1.D#1998.year 1.D#1999.year 1.D#2000.year 1.D#2001.year 1.D#2002.year 1.D#2003.year 1.D#2004.year /// 1.D#2005.year ///
	 1.D#2006.year 1.D#2007.year 1.D#2008.year 1.D#2009.year 1.D#2010.year 1.D#2011.year 1.D#2012.year 1.D#2013.year ///
	 c.X#i.year c.DX#i.year if inrange(tot_p01,150,450) & pop_mismatch20==0 & lights_diff<20 [aw=aw150], a(pca01_id year#stdt) vce(cluster stdtbk) 
	 
cap drop y beta lci uci
gen y = _n+1997
replace y = . if y>2013
gen beta = .
gen lci = .
gen uci = .

forvalues y = 1998/2013 {
	cap replace beta = _b[1.D#`y'.year] if y==`y'
	cap replace lci =  _b[1.D#`y'.year] - 1.96*_se[1.D#`y'.year] if y==`y'
	cap replace uci =  _b[1.D#`y'.year] + 1.96*_se[1.D#`y'.year] if y==`y'
}
cap replace beta = 0 if y==2005

twoway (rspike uci lci y, lstyle(ci) lw(medium) lcolor(gs10)) ///
	(scatter beta y, msize(medsmall) mc(midblue) msymbol(O)) ///
	, ///
	yline(0, lcolor(black) lw(thin)) /// 
	ytitle("Year-specific RD coefficient (95% CI)", size(medlarge)) xtitle("Nighttime brightness year", size(medlarge)) ///
	title("District-year FEs, 150-person bandwidth", color(black) size(large)) /// 
	graphregion(lcolor(white)) graphregion(color(white)) plotregion(fcolor(white)) ///
	ylabel(-0.1 0 0.1 0.2 0.3 0.4,nogrid angle(0) labsize(medlarge)) yscale(range(-0.1 0.43)) ///
	xlabel(1998(1)2013, labsize(medlarge) angle(45)) ///
	 legend(off) //aspect(0.4)

graph export "$results/dd figs/diff_in_disc_lights_dt_150.pdf", replace 
	 
}

********************************************************************************
********************************************************************************
	
** 4. Difference-in-Discontinuities, heterogeneous RD in hours/day of power supply
{
use "${panel}/panel_dataset_dd.dta", clear
keep if vplan4<11 & corr_state==1 & sample==1 & sample_h==1 & sing_h==1 & pop_non_zero==1 
drop pct_06-vd_ame_d_spcl
drop vd_com_d_newspap-vd_ame_d_cv_hall
drop work_main_ag_p-pct_irr_twell pct_irr_irr_twell
drop p_300_X-band300_p_200
drop bin25_p-within_border_100k
drop vplan vplan2 vplan3 vdpr vdpr2 vdpr3
drop treat_x_post dd_* stateyear
drop lights_mean_hat lights_max_hat
drop no_hh tot_p tot_m tot_f vd01_id conc_id v_id_hpca11 vd11_id c_code01 v_id_master ///
	h?_id h?v_id f_area totp_h? count_h? maxp_h? grad_max grad_mean
merge 1:1 pca01_id st_code dt_code vplan4 vdpr4 year using ///
	"${panel}/panel_dataset_lights_allyears.dta"
egen max_merge = max(_merge), by(pca01_id)
keep if max_merge==3
tab year _merge
assert _merge==3 if year==2001 | year==2011
drop _merge max_merge

order pca01_id year state st_code district dt_code stdt vplan4 vdpr4 lights* vd_pwr*
preserve
keep if inlist(year,2001,2011)
keep pca01_id bk_code-popbin_2700_max
duplicates drop
unique pca01_id
assert r(unique)==r(N)
tempfile vars_constant
save `vars_constant'
restore
drop bk_code-popbin_2700_max
merge m:1 pca01_id using `vars_constant'
assert _merge==3
drop _merge
sort pca01_id year

preserve
use "${panel}/panel_dataset_full.dta", clear
keep pca01_id lights_max????_hat
rename lights_max????_hat lights_max_hat????
reshape long lights_max_hat, i(pca01_id) j(year)
tempfile lights_max_hat
save `lights_max_hat'
restore
merge 1:1 pca01_id year using `lights_max_hat', keep(1 3)
drop if tot_p01>1000

egen temp1 = max(lights_max*(year==2001)), by(pca01_id)
egen temp2 = max(lights_max*(year==2011)), by(pca01_id)
gen lights_diff = abs(temp2-temp1)
assert lights_diff!=.

egen stdtbk = group(stdt bk_code)
gen X = (tot_p01-300)
assert X!=.
gen D = X>=0
gen DX = D*X

gen aw150 = abs(150-abs(300-tot_p01))
gen aw100 = abs(100-abs(300-tot_p01))


preserve
use "$panel/panel_dataset_full.dta", clear
egen temp1 = mean(vdp_pwr_h_all_avg) if vdp_pwr_h_all_avg_11>0 & vdp_pwr_h_all_avg_11!=., by(st_code dt_code)
egen HRS_all_wide = mode(temp1), by(st_code dt_code)
egen temp2 = mean(vdp_pwr_h_dom_avg) if vdp_pwr_h_dom_avg_11>0 & vdp_pwr_h_dom_avg_11!=., by(st_code dt_code)
egen HRS_dom_wide = mode(temp2), by(st_code dt_code)
gen HRS_all = HRS_all_wide>=10 & HRS_all_wide!=.
gen HRS_dom = HRS_dom_wide>=10 & HRS_dom_wide!=.
replace HRS_all = . if HRS_all_wide==.
replace HRS_dom = . if HRS_dom_wide==.
keep st_code dt_code HRS_all HRS_dom
duplicates drop
tempfile hrs
save `hrs'
restore
merge m:1 st_code dt_code using `hrs', nogen

cap drop H
gen H = HRS_all

unique st_code dt_code if H==1 & inrange(tot_p01,200,400) & pop_mismatch20==0 & lights_diff<20
unique st_code dt_code if H==0 & inrange(tot_p01,200,400) & pop_mismatch20==0 & lights_diff<20
unique st_code if H==1 & inrange(tot_p01,200,400) & pop_mismatch20==0 & lights_diff<20
unique st_code if H==0 & inrange(tot_p01,200,400) & pop_mismatch20==0 & lights_diff<20
tab st_code if H==1 & inrange(tot_p01,200,400) & pop_mismatch20==0 & lights_diff<20
tab st_code if H==0 & inrange(tot_p01,200,400) & pop_mismatch20==0 & lights_diff<20

 // District FEs, 100 bw
reghdfe lights_max ///
	 1.D#1998.year#1.H 1.D#1999.year#1.H 1.D#2000.year#1.H 1.D#2001.year#1.H 1.D#2002.year#1.H 1.D#2003.year#1.H 1.D#2004.year#1.H /// 1.D#2005.year ///
	 1.D#2006.year#1.H 1.D#2007.year#1.H 1.D#2008.year#1.H 1.D#2009.year#1.H 1.D#2010.year#1.H 1.D#2011.year#1.H 1.D#2012.year#1.H 1.D#2013.year#1.H ///
	 1.D#1998.year#0.H 1.D#1999.year#0.H 1.D#2000.year#0.H 1.D#2001.year#0.H 1.D#2002.year#0.H 1.D#2003.year#0.H 1.D#2004.year#0.H /// 1.D#2005.year ///
	 1.D#2006.year#0.H 1.D#2007.year#0.H 1.D#2008.year#0.H 1.D#2009.year#0.H 1.D#2010.year#0.H 1.D#2011.year#0.H 1.D#2012.year#0.H 1.D#2013.year#0.H ///
	 c.X#i.year c.DX#i.year if inrange(tot_p01,200,400) & pop_mismatch20==0 & lights_diff<20 [aw=aw100], a(pca01_id year#stdt) vce(cluster stdtbk) 
	 
cap drop y beta_hi lci_hi uci_hi beta_lo lci_lo uci_lo
gen y = _n+1997
replace y = . if y>2013
gen beta_hi = .
gen lci_hi = .
gen uci_hi = .
gen beta_lo = .
gen lci_lo = .
gen uci_lo = .

forvalues y = 1998/2013 {
	cap replace beta_hi = _b[1.D#`y'.year#1.H] if y==`y'
	cap replace lci_hi =  _b[1.D#`y'.year#1.H] - 1.96*_se[1.D#`y'.year#1.H] if y==`y'
	cap replace uci_hi =  _b[1.D#`y'.year#1.H] + 1.96*_se[1.D#`y'.year#1.H] if y==`y'
	cap replace beta_lo = _b[1.D#`y'.year#0.H] if y==`y'
	cap replace lci_lo =  _b[1.D#`y'.year#0.H] - 1.96*_se[1.D#`y'.year#0.H] if y==`y'
	cap replace uci_lo =  _b[1.D#`y'.year#0.H] + 1.96*_se[1.D#`y'.year#0.H] if y==`y'
}
cap replace beta_hi = 0 if y==2005
cap replace beta_lo = 0 if y==2005

cap drop y_hi y_lo
gen y_hi = y-0.17
gen y_lo = y+0.17

twoway ///
	(rspike uci_hi lci_hi y_hi, lstyle(ci) lw(medium) lcolor(gs10)) ///
	(rspike uci_lo lci_lo y_lo, lstyle(ci) lw(medium) lcolor(gs10)) ///
	(scatter beta_hi y_hi, msize(medium) mc(midblue) msymbol(O)) ///	
	(scatter beta_lo y_lo, msize(medium) mc(maroon) mfc(white) msymbol(O)) ///
	, ///
	yline(0, lcolor(black) lw(thin)) /// 
	ytitle("Year-specific RD coefficient (95% CI)", size(medlarge)) xtitle("Nighttime brightness year", size(medlarge)) ///
	title("District-year FEs, 100-person bandwidth", color(black) size(large)) /// 
	graphregion(lcolor(white)) graphregion(color(white)) plotregion(fcolor(white)) ///
	ylabel(-0.1(0.1)0.5,nogrid angle(0) labsize(medlarge)) yscale(range(-0.1 0.53)) ///
	xlabel(1998(1)2013, labsize(medlarge) angle(45)) ///
	legend(order(3 "{&ge}10 hours/day of power   " ///
	             4 "<10 hours/day of power" ) ///
			c(2) pos(6) size(medlarge))

graph export "$results/dd figs/diff_in_disc_lights_dt_100_hours.pdf", replace 
	 

 // District FEs, 150 bw
reghdfe lights_max ///
	 1.D#1998.year#1.H 1.D#1999.year#1.H 1.D#2000.year#1.H 1.D#2001.year#1.H 1.D#2002.year#1.H 1.D#2003.year#1.H 1.D#2004.year#1.H /// 1.D#2005.year ///
	 1.D#2006.year#1.H 1.D#2007.year#1.H 1.D#2008.year#1.H 1.D#2009.year#1.H 1.D#2010.year#1.H 1.D#2011.year#1.H 1.D#2012.year#1.H 1.D#2013.year#1.H ///
	 1.D#1998.year#0.H 1.D#1999.year#0.H 1.D#2000.year#0.H 1.D#2001.year#0.H 1.D#2002.year#0.H 1.D#2003.year#0.H 1.D#2004.year#0.H /// 1.D#2005.year ///
	 1.D#2006.year#0.H 1.D#2007.year#0.H 1.D#2008.year#0.H 1.D#2009.year#0.H 1.D#2010.year#0.H 1.D#2011.year#0.H 1.D#2012.year#0.H 1.D#2013.year#0.H ///
	 c.X#i.year c.DX#i.year if inrange(tot_p01,150,450) & pop_mismatch20==0 & lights_diff<20 [aw=aw150], a(pca01_id year#stdt) vce(cluster stdtbk) 
	 
cap drop y beta_hi lci_hi uci_hi beta_lo lci_lo uci_lo
gen y = _n+1997
replace y = . if y>2013
gen beta_hi = .
gen lci_hi = .
gen uci_hi = .
gen beta_lo = .
gen lci_lo = .
gen uci_lo = .

forvalues y = 1998/2013 {
	cap replace beta_hi = _b[1.D#`y'.year#1.H] if y==`y'
	cap replace lci_hi =  _b[1.D#`y'.year#1.H] - 1.96*_se[1.D#`y'.year#1.H] if y==`y'
	cap replace uci_hi =  _b[1.D#`y'.year#1.H] + 1.96*_se[1.D#`y'.year#1.H] if y==`y'
	cap replace beta_lo = _b[1.D#`y'.year#0.H] if y==`y'
	cap replace lci_lo =  _b[1.D#`y'.year#0.H] - 1.96*_se[1.D#`y'.year#0.H] if y==`y'
	cap replace uci_lo =  _b[1.D#`y'.year#0.H] + 1.96*_se[1.D#`y'.year#0.H] if y==`y'
}
cap replace beta_hi = 0 if y==2005
cap replace beta_lo = 0 if y==2005

cap drop y_hi y_lo
gen y_hi = y-0.17
gen y_lo = y+0.17

twoway ///
	(rspike uci_hi lci_hi y_hi, lstyle(ci) lw(medium) lcolor(gs10)) ///
	(rspike uci_lo lci_lo y_lo, lstyle(ci) lw(medium) lcolor(gs10)) ///
	(scatter beta_hi y_hi, msize(medium) mc(midblue) msymbol(O)) ///	
	(scatter beta_lo y_lo, msize(medium) mc(maroon) mfc(white) msymbol(O)) ///
	, ///
	yline(0, lcolor(black) lw(thin)) /// 
	ytitle("Year-specific RD coefficient (95% CI)", size(medlarge)) xtitle("Nighttime brightness year", size(medlarge)) ///
	title("District-year FEs, 150-person bandwidth", color(black) size(large)) /// 
	graphregion(lcolor(white)) graphregion(color(white)) plotregion(fcolor(white)) ///
	ylabel(-0.1(0.1)0.5,nogrid angle(0) labsize(medlarge)) yscale(range(-0.1 0.53)) ///
	xlabel(1998(1)2013, labsize(medlarge) angle(45)) ///
	legend(order(3 "{&ge}10 hours/day of power   " ///
	             4 "<10 hours/day of power" ) ///
			c(2) pos(6) size(medlarge))

graph export "$results/dd figs/diff_in_disc_lights_dt_150_hours.pdf", replace 
	 
}

********************************************************************************
********************************************************************************

